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I.  INTRODUCTION 


Although  many  computer  programs  for  the  solution  of  space-charge 
electron  flow  problems  exist j1-1*  there  is  a  need  for  a  program  which  is 
user  friendly,  well  documented,  portable,  versatile,  and,  most  impor¬ 
tantly,  available  at  no  cost  to  the  scientific  and  industrial  community. 

This  report  describes  the  project  to  design  and  implement  a 
deformable  triangle  mesh  computer  program  which  solves  the  space-charge 
electron  flow  problem.  The  computer  program  was  implemented  on  an  HP- 
1000  computer  at  Teledyne  MEC  in  Palo  Alto,  California. 


The  analysis  of  the  space-charge  flow  problem  consists  of  two 


parts: 

1.  Solving  Poisson's  equation 

2.  Solving  the  equation  of  motion  for  electrons  derived  from  the 
Lorentz  force  equation 

To  solve  Poisson's  equation,  the  set  of  codes  developed  by  Halbach  and 
Holsinger,5  "Poisson  Group  Programs,”  was  used.  These  codes  were  chosen 
because  of  their  reliable  results,  ability  to  solve  both  electrostatic 
and  magnetostatic  problems,  and  the  option  to  change  the  permittivity  or 
permeability  in  different  regions.  The  solution  of  the  Lorentz  force 
equation  follows  work  published  by  R.  True,6  and  by  Caplan  and  Thoring- 
ton,7  with  some  modifications. 

For  this  project,  the  following  objectives  have  been  achieved: 

1.  The  computer  program  is  user  friendly.  In  most  deformable 
mesh  programs,  the  generation  of  the  analysis  mesh  is 
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particularly  troublesome.  To  generate  the  triangular  mesh,  a 
nodal  diagram  must  be  constructed,  and  the  nodal  indices  of 
all  boundary  nodes  must  be  determined.  The  nodal  indices,  as 
well  as  the  corresponding  physical  boundary  points,  must  be 
supplied  as  input  to  the  program.  This  is  a  somewhat  time- 
consuming  process  which  requires  some  experience.  In  order  to 
avoid  this  problem,  the  Poisson  group  programs  include  an 
automatic  mesh  generation  program  called  AUTOMESH.  AUTOMESH 
constructs  the  "logical*1  mesh  and  generates  (x,  y)  coordinate 
data  which  describes  the  physical  boundary. 

The  program  is  well  documented  and  portable.  The  computer 
code  is  written  in  ANSI  Standard  FORTRAN  77,  so  the  program 
may  be  transferred  between  computers.  Since  the  program  is 
well  documented,  it  is  more  readable  and,  hence,  easier  to 
understand,  use,  and  modify. 

The  computer  code  is  very  versatile.  For  two-dimensional 
problems  In  rectangular  or  cylindrical  geometry,  it  will 
solve: 

a.  Space-charge  flow  problems 

b.  Electrostatic  problems 

c.  Nonlinear  magnetostatic  problems  with  saturated  permea¬ 
bilities. 

In  the  future,  the  ability  to  solve  for  the  electrostatic  and 
nonlinear  magnetostatic  fields  from  the  same  mesh  will  be 
Incorporated  into  the  code. 
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II.  SPACE-CHARGE  FLOW 


2.1  The  Space-Charge  Flow  Problem 


Two  aets  of  basic  equations  are  required  in  space-charge  flow. 

The  first  set  consists  of  the  equations  obeyed  by  the  static  fields, 

assuming  that  the  charge  and  current  densities  are  known.  The  second 

set  predicts  the  motion  of  the  particles  in  given  electrostatic  and 

magnetostatic  fields.  The  first  set  can  be  deduced  from  Maxwell's 

equations.  The  second  set  can  be  deduced  from  the  Lorentz  force  law, 

with  the  help  of  the  field  equations. 

♦  + 

The  electric  field  E  and  the  magnetic  field  B  must  obey  Maxwell's 
equations,  which  take  the  form,  for  steady-state  flows  in  a  vacuum, 


V  •  (el)  -  f 


7  x  £  -  0 


V  •  B  *  0 


7  x  B  -  -y  i 
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where  p  and  i  are  charge  and  current  densities. 

From  Eq.  1,  it  is  possible  to  define  an  electric  scalar  potential 
V  so  that 


E  «  -  7V 
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With  the  V  of  Eq.  5,  it  is  possible  to  derive  Poisson's  equation  from 
Eq.  1. 

7  .  (erW)  -  -  (6) 

o 

The  second  set  of  equations  which  govern  the  motion  of  the  particles  in 
the  field  are  deduced  from  the  Lorentz  force  law.  The  force  F  on  a 
particle  in  an  electromagnetic  field  is 

F  -  q(E  +  v  x  b)  (7) 

♦ 

where  v  is  the  velocity  of  the  particle.  From  Eq.  7,  the  acceleration 
and  velocity  of  an  electron  are  determined  in  terms  of  the  fields  in  the 
system.  The  fields  which  affect  the  electron  motion  consist  of  two 
parts: 

1.  Those  set  up  by  the  electrodes  and  external  magnetic  sources 

2.  Those  due  to  the  electrons  in  the  beam. 

Although  Eq.  4  is  the  exact  equation  obeyed  by  the  magnetic  field,  for 
nonrelativistlc  velocities,  it  is  usually  reasonable  to  neglect  the 
self-magnetic  field  due  to  the  beam  itself.  The  self  magnetic  forces  on 
a  particle  are  usually  small  compared  with  the  space-charge  forces.8 


Solution  of  the  Space-Charge  Flow  Problem 


A  deformable  mesh  electron  gun  analysis  program  solves  the  space- 
charge  flow  problem  by  computing  in  a  self-consistent  way  the  poten¬ 
tials,  current  density,  electron  trajectories,  and  space-charge  distri¬ 
bution. 

The  general  method  of  solution  is  shown  in  the  flow  chart  of  Fig. 

1.  The  iterative  process  can  be  described  as  follows: 

1.  Solve  for  the  potential  at  each  mesh  point  given  the  boundary 
conditions,  and  space-charge  density  at  each  mesh  point  by 
using  Poisson' 8  equation. 


V  •  (erW)  -  -  (8) 

o 

2.  Calculate  the  current  density  along  the  cathode  and  the 
emitted  current  using  Child's  law. 


3.  Compute  the  electron  trajectories  from  the  Lorentz  force 
equation, 

F  -  q(E  ♦  v  x  b)  (10) 


given  the  electric  and  magnetic  fields. 


A.  Assign  charge  to  each  node  the  electron  passes  near. 

5.  Repeat  steps  1  through  A  until  the  process  is  self-consistent. 


INPUT  DATA 


III.  SOLUTION  OF  POISSON'S  EQUATION 


3.1  Choice  of  Method 


The  nonlinear,  two-dimensional  Poisson  equation  is  solved  by  the 
finite  difference  method  using  a  nonuniform  triangular  mesh.  The  advan¬ 
tages  of  the  triangular  mesh  for  the  numerical  solution  of  partial 
differential  equations  in  two  dimensions  were  first  pointed  out  in  1943 
by  Courant.9  Later  papers  by  True10  and  Cattellno11  give  examples 
demonstrating  the  advantages  and  disadvantages  of  this  technique.  The 
most  important  advantage  of  the  deformable  triangle  mesh  compared  to  the 
rectangular  mesh  is  that  the  mesh  density  can  be  selectively  varied 
within  the  analysis  region.  The  mesh  density  can  be  increased  in 
regions  where  the  mesh  spacing  must  be  small  compared  to  the  features 
being  modeled,  and  the  density  can  be  decreased  in  regions  where  the 
field  gradients  sre  small.  Therefore,  there  is  no  need  to  increase  the 
number  of  nodes  and,  hence,  the  number  of  calculations  remain  the  same. 

Other  advantages  of  triangular  mesh  analysis  over  the  rectangular 
mesh  analysis  include: 

1.  Boundary  positions  are  fixed  and  mesh  nodes  are  positioned  to 
lie  exactly  on  the  boundaries.  This  is  not  the  case  in  the 
rectangle  mesh  analysis. 

2.  Errors  in  the  potential  are  generally  reduced  because  each 
mesh  point  has  six  neighbors  instead  of  four. 

3.  Since  a  high  density  "cathode  mesh"  is  generated  near  the 
cathode,  a  large  improvement  in  the  accuracy  of  calculated 
values  for  emission  current  density  and  perveance  is  obtained. 
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As  sight  be  expected,  there  are  also  some  disadvantages  in  the 
triangular  mesh  analysis  technique.  The  use  of  a  triangular  mesh 
results  in  considerably  more  complex  programming  logic  as  compared  with 


the  rectangular  analysis.  Other  disadvantages  include: 

1.  More  computer  storage  is  needed  per  mesh  point.  This  is  due 
to  saving  three  coupling  coefficients  and  (x,  y)  coordinates 
of  each  mesh  point  in  addition  to  the  voltage  and  space-charge 
density  at  that  point. 

2.  Increase  in  time  needed  to  compute  the  trajectories  caused  by 
the  need  to  search  for  the  nearest  mesh  point  after  every  time 
step. 

Thus,  the  increased  accuracy  obtained  is  at  the  price  of  producing 
a  more  complicated  program.  However,  these  effects  are  not  drastic  and, 
in  cases  such  as  the  gridded  gun  where  accuracy  is  needed  in  one  region 
and  unnecessary  in  other  regions,  the  triangular  mesh  analysis  is  the 
best  way  to  proceed. 

3.2  Description  of  the  Method 

The  detailed  derivation  of  the  equations  for  the  finite  difference 
approximations  on  a  general  triangle  mesh  are  found  in  Winslow.12  The 
equations  that  are  used  in  the  computer  program  are  summarized  in  Sec¬ 
tions  3.3  and  3.4. 

Poisson's  equation. 
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la  to  ba  aolvad  ovar  a  region  R,  where  the  permittivity  is  a  positive 
function  of  the  rectangular  coordinates  (x,  y),  and  the  charge  density  p 
la  a  given  function  of  (x,  y). 

The  baaic  aaaunptlona  of  the  finite  difference  method  are: 

1.  The  boundaries  and  Interfaces  of  the  region  R  are  approximated 
by  straight-line  segments. 

2.  The  region  is  triangulated. 

3.  The  values  of  the  potential  V  are  defined  at  triangle  vertices 
and  V  la  assumed  to  vary  linearly  over  each  triangle. 

4.  The  charge  density  and  permittivity  are  assumed  to  be  constant 
over  each  triangle. 


3.3  Difference  Equations 


The  triangulation  used  in  this  method  is  topologically  regular; 
six  triangles  swet  at  every  interior  mesh  point.  Consider  an  interior 
mesh  point  in  the  triangular  mesh,  as  shown  in  Fig.  2.  A  secondary  mesh 
of  12  aides  whose  vertices  are  alternately  the  centroids  of  the  six 
adjacent  triangles  and  the  midpoints  of  the  six  adjacent  sides  is 
defined  with  respect  to  the  primary  triangle  mesh.  This  is  the  shaded 
region,  as  shown  in  Fig.  2.  The  secondary  mesh  element  is  comprised  of 
one-third  of  the  area  of  each  of  the  six  primary  mesh  triangles  sharing 
that  vertex,  so  that  each  triangle  of  area  A  is  divided  into  three  equal 
quadrilaterals  of  area  a  *  A/3.  This  quadrilateral  area,  a,  will  be 
used  to  compute  the  source  terms  S  of  Eq.  12. 
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n? 

m 


RES 


For  any  internal  node,  the  finite  difference  analog  of  Eq.  11  la 


I  -  V0)  +  So  -  0 


where  the  subscript  o  represents  the  point  at  which  the  calculation  is 
made,  and  1  represents  the  six  neighboring  vertices.  The  tera  w^  is  the 
coupling  coefficient  for  the  line  joining  the  vertex  i  and  the  center, 
and  depends  only  on  the  nature  of  the  two  triangles,  having  this  as  a 
coBBon  side.  Froa  Winslow,13  the  coupling  coefficient  is  given  as 


'i  •  T  (s 


1+1/2  +  ,  i-1/2 

cot  0  ♦  e  cot  0 

r  r 


where  the  angles  0  and  0  ate  defined  in  Fig.  2. 


Note  that  the  coupling  between  two  points  (x^,  yj)  and  (x^,  y^}  is 


syametric,  so  that 


w12  -  w2J 


Since  the  evaluation  of  trlgonoaetrlc  functions  with  a  computer  is 
inefficient,  an  alternate  expression  for  coupling  coefficients  Is 

derived  in  True14  and  is  given  as 

* 

1  1+1/2  (Vl  ~  Xo^Xi+l  Xi^  *  ^yi+l  ~  yo^yi+l  ~  yj) 

” 2 1  e,r  L  ^yi '  yo^xi+i "  xoJ  +  [xo _  Xi^yi+i " 


i-1/2  pXi-l  ~  XoJ(xi-l  "  *  (yi-l  ~  yo)(yi-l  1  yl] 

r  lyi_i  -  y0J(xi  ~  X0J  +  tx0  ~  xi-i^yi  " 
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The  source  ttri  S  of  Eq*  12  is  conputed  ss 


S0  “  Si+l/2  *i+l/2 


where  Sj+j^2  ar*d  ai+i/2  are  tlie  8°urce  density  and  quadrilateral  area, 
respectively,  for  triangle  1  +  1/2  of  Fig.  2.  The  source  density  of 
triangle  i  +  1/2  is  defined  by 


i+1/2 


i+1/2 


and  Eq.  16  becomes 


S  -  V  °l+l/2  *1+1/2 
°  i-1  6o 


where  Pi+1/2  is  the  char®e  density  of  triangle  i  +  1/2. 

For  each  mesh  point,  Eq.  12  is  solved  by  the  iterative  method  of 
successive  overrelaxation.15  The  difference  formula,  which  must  be 
satisfied  at  each  interior  point,  is  given  by 


?  Pi+l/2  *1+1/2 


l  V  w  +  |  l 
1-1  1-1 
6 

l 

i-1 


Introducing  the  overrel axat ion  parameter  RH0A1R  (0  <  RHOAIR  <  2),  we 


V./V.  •*. 


where  RHOAIR  is  a  program  parameter  which  is  used  to  accelerate  the 
convergence  for  systems  that  are  convergent  by  the  Gauss-Seidel  tech¬ 
nique.  For  each  iteration  cycle  of  Eq.  20,  the  mesh  points  are  swept 
in  sequence  in  which  the  nearest-neighbor  values  v",n+*  represent 
if  it  already  has  been  calculated,  or  V°  if  it  has  not. 

The  program  continues  to  iterate  until 


ABS(v"+1  -  V®)  <  EPSILA 


where  EPSILA  is  the  convergence  criterion  for  the  potential  solution, 
which  is  typically  set  to  10-6.  EPSILA  is  also  a  program  parameter. 

3.4  Extension  to  Problems  in  the  Z-R  Plane 

Assuming  the  permittivity  e  is  not  constant  throughout  the  region, 
Poisson '8  equation  may  be  expressed  as 


y  .  (erW)  -  f 
o 


Evaluating  Eq.  22  in  cylindrical  coordinates  and  assuming  no  variation 
in  the  8  direction, 


1  8  f  3v1  A  1  3  T  3V 1  -p 

r  3r  j^rCr  3rJ  r  3z  |^rer  dzj  " 


.-.V.V 
,•  V  V  ; 

» 

V  V  v 

V  J 
■>>; 
'v‘v  V 

o.Vj 


>>V: 

8S® 


>'.V 

a.V 

,V/% 


cv*-'V.v /y./. s*  1 


>  ■  ^  *  _>  Va  ■ 


Now,  considering  the  radial  weighting  of  the  charge  density  in  Eq.  26b, 
the  source  tern  SQ  of  Eq.  16  is  modified  for  cylindrically  symmetric 
problems  as 


S 

o 


(ri+l/2  Pl+l/2^ 

e  ai+l/2 

o 


(29) 


where  ri+i/2  *8  t*'e  overage  radius  of  a  quadrilateral  at  vertex  rQ  and 
is  given  by 


_7  +  5 

ri+l/2  12  r  12  ro 


(30) 


Therefore,  it  is  possible  to  use  virtually  the  same  programming  logic 
used  in  the  x  -  y  plane  problem  by  interchanging  "z"  as  "x"  and  “r”  as 
"y"  and  by  modifying  the  coupling  coefficients  and  the  source  terms,  as 
described  by  Eq.  27  and  29,  respectively. 


IV..  CURRENT  EMISSION 


Once  Poisson's  equation  has  been  solved,  the  current  emitted  from 
each  cathode  point  is  calculated  by  the  Child-Langmuir  space-charge  law 
for  the  planar  diode. 17,18  After  the  current  has  been  determined, 
charge  is  assigned  in  the  near  cathode  region  and  the  perveance  is 
calculated. 

4.1  Planar  Diode  Approximation  Near  the  Emitting  Surface 

The  emitting  surface  consists  of  a  set  of  elementary  planar 
diodes.  Current  emission  is  calculated  from  each  elementary  diode  by 
Childe's  law  based  on  the  voltages  at  nodes  in  front  of  the  cathode 
surface. 

Near  the  cathode  surface,  an  approximately  rectangular  mesh  is 
generated,  as  shown  in  Fig.  3.  The  spacing  "z"  of  the  rectangular  mesh 
at  node  m  is  determined  by 

z  -  (a  +  b)/(FRT  *  NAC)  (31) 
m 

where  z,  a,  and  b,  are  defined  in  Fig.  3.  The  program  parameter  FRT 
controls  the  nearness  of  the  nodes,  and  the  program  parameter  NAC  is  the 
number  of  nodes  in  front  of  the  emitting  surface. 

It  has  been  determined  empirically  that,  for  curved  emitting 
surface  being  approximated  by  a  set  of  planar  diodes,  two  criteria 
should  simultaneously  be  satisifed:19 
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1.  The  change  In  the  angle  of  the  line  of  nodes  in  the  rectangu¬ 
lar  mesh  from  cathode  point  to  cathode  point  should  be  less 
than  4  degrees. 

2.  The  program  parameter  FRT,  which  controls  the  nearness  of  the 
nodes,  should  be  greater  than  8. 

The  user  must  be  careful  that  these  two  conditions  are  satisfied,  since 
the  potentials  deviate  substantially  from  the  planar  diode  theory  when 
they  are  not  close  to  the  emitting  surface* 

4.2  Calculation  of  Emission 

In  a  space-charge  limited  parallel  plane  diode,  the  potential 
varies  as  the  four-thirds  power  of  distance  from  the  emitting  surface, 

V(x)  -  Cx4/3  (32) 

where  C  is  a  constant  and  x  is  the  distance  normal  to  the  surface.  The 
current  density  for  a  planar  diode  is  given  as 

3/2 

J  -  - - kC*'*  (33) 

x 

where 

k  -  4/9  e  /2 rT 
o 

is  the  permittivity  of  free  space 
n  is  the  electron  charge  to  mass  ratio 

Near  the  cathode  surface  is  the  previously  discussed  rectangular 
mesh,  which  is  NAC  rectangles  thick,  as  shown  in  Fig.  3.  At  each  node 


on  the  emitting  surface,  the  constant  C  Is  determined  from  Eq.  32  for 
each  of  the  HAC-1  nodes  in  front  of  the  cathode.  The  average  C  of  these 
RAC-1  nodes  is  then  taken  as  the  C  for  that  emitting  node.  Letting  m 


l 

\ 


\ 


denote  nodes  along  the  cathode  and  t  denote  nodes  normal  to  the  cathode, 
the  diode  constants  Cm  are  computed  from 


C»- 


NAC-1 

I 

.  i-1 


(Ml 


473 


(NAC-1) 


(34) 


At  this  point,  the  program  laterally  averages  the  calculated  C's  with 
its  neighbors  to  suppress  oscillation  in  C  along  the  emitting  surface. 
This  oscillation  arises  from  nonuniform  emission  from  the  cathode  sur¬ 
face.  Let  C„  be  the  constant  for  the  emitting  node  currently  being 
averaged,  be  the  constant  for  the  emitting  node  below,  and  be 
the  constant  for  the  emitting  node  above.  C  is  calculated  by  laterally 
averaging  as  follows: 


C 

m 


(C  .  +  2C  +  C  ..) 
'  m-1  m  m+1 


(35) 


If  m  is  an  axial  node,  the  lateral  averaging  is  performed  as 


(Vi + 


(36) 


Por  upper  edge  nodes,  Cm  is 


C 

m 


+Cm-l) 
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(37) 


In  the  first  few  iterations,  C  will  be  too  large,  since  space-charge 

m 

depression  has  not  yet  occurred  in  front  of  the  emitting  surface.  To 
correct  this  problem,  the  program  next  averages  the  diode  constant  C 

m 

with  the  C°ld  from  the  previous  cycle  by  the  equation 
m 


c  -  cold  +  y(c  -  cold) 

m  m  v  m  m  J 


(38) 


where  Y  (or  CGAMHA)  is  a  program  parameter.  A  typical  value  for  CGAMMA 
should  be  less  than  0.3,  and  it  should  be  made  smaller  for  a  thinner 
rectangular  mesh  near  the  cathode.  Although  a  small  CGAMMA  will  prevent 
unstable  numerical  oscillations  in  the  current  between  cycles,  choosing 
a  CGAMMA  too  small  will  make  the  program  go  through  too  many  cycles 
before  the  solution  converges. 

If  negative  potentials  exist  in  front  of  the  emitting  surface,  C^ 
will  become  negative.  To  avoid  this  situation,  the  program  will  compute 


C 

m 


-old 

SSSF  •  C 

m 


(39) 


where  SSSF  is  a  program  parameter  which  decreases  the  emitted  cur¬ 
rent,20  and  is  typically  set  to  0.2. 

The  temperature  limited  value  of  C  is  determined  from  Eq.  33  as 


C 


temperature  limited 


^FJLIMI 


) 


2/3 


(40) 


and  is  given  by  the  user 


4 


S 


±1 


where  FJLIMI  is  the  current  density  in  — 

_  m 

If  C  is  greater  than  this  value,  then  C  is  recomputed  as 
Hi  is 


C  m  C 

m  temperature  limited 


where  the  temperature  limited  value  of  the  current  density  is  given  to 
the  program  by  the  user. 


4.2.1  Temperature  Limited  Emission 


The  current  emitted  from  the  cathode  may  be  set  to  a  fixed 
value.  The  program  parameter,  CURREN,  is  the  amount  of  current  in  amps 
emitted  by  the  cathode.  If  CURREN  is  greater  than  zero,  the  program 
will  compute  the  current  density,  JCUR,  at  the  emitting  mode  m  by 


CURREN  .... 

JCURm  “  NODES  CAREA  (  ) 

m 

where  CAREAm  is  given  in  Eq.  53  and  NODES  is  the  number  of  emitting 
points.  Once  JCUR^  has  been  computed  at  the  emitting  node,  the  emission 
constant,  C  ,  is  recomputed  from  Eq.  33  as 

IS 


!.-(=■)' 


and  the  program  proceeds  normally. 


4.3  Program  Convergence 


The  computer  program  decides  convergence  in  the  iterative  pro¬ 
cedure  by  the  relative  change  in  between  major  current  cycles  (see 
Fig.  1). 

The  problem  is  assumed  to  have  converged  when 


DC  BAR  <  ECBAR 


where 


DC  BAR 


(44) 


(45) 


and  ECBAR  is  a  program  parameter  specified  by  the  user. 

The  user  may  also  limit  the  number  of  major  current  cycles  to  be 
performed  by  specifying  the  program  parameter  MAXCYC,  which  is  the 
maximum  number  of  current  cycles  to  calculate. 


4.4  Charge  Assignment  in  the  Wear-Cathode  Region 

Near  the  cathode,  it  is  assumed  that  planar  diode  theory 
applies.  Charge  assignment  to  nodes  1  through  NAC-1  of  the  rectangular 
mesh  adjacent  to  the  cathode  are  obtained  by  multiplying  the  charge 
density  calculated  at  the  node  by  the  cell  area  times  a  special  weight¬ 
ing  function.  This  technique  of  providing  high  accuracy  near  the  cath¬ 
ode  leads  to  a  significant  improvement  in  current  calculations. 

For  a  cathode  mesh  being  approximately  rectangular,  one  may  assume 
that  the  problem  is  one-dimensional.21  Poisson's  equation  for  the  one- 
dimensional  case  reduces  to 
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where  x  is  the  direction  normal  to  the  emitting  surface.  The  space- 
charge  density  for  s  planar  diode  can  be  obtained  from  Eq.  46  and  i6 
given  as 


p(x) 


(47) 


Using  Eq.  47,  the  source  terms  for  nodes  1  through  NAC-1  in  the  near 
cathode  region  are  computed  as 


S 


t 


-M  p(x)  ADOD£ 


(48) 


where  is  a  special  weighting  function  and  is  given  by 


M 


l 


2 

4 


(1  ♦  1)4/3  -  2i4/3 


,-2/3 


(49) 


ADOD^  is  the  dodecagon  area  of  the  1th  node  adjacent  to  the  emitting 
point  and  is  given  as 


6 

ADOD  l  a.  (50) 

i-1  1 

where  a^  is  the  quadrilateral  area  for  triangle  i  of  Fig.  2. 

According  to  True,22  the  use  of  the  special  weighting  function 
is  introduced  because  the  potentials  calculated  by  the  finite-difference 
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method  do  aot  agro*  with  the  theoretical  values  predicted  by  Child's 
law.  This  discrepancy  Is  partly  due  to  the  truncation  error  in  the 
finite  difference  approximation  and  that  the  charge  density  is  assumed 
aot  to  vary  within  the  cell. 


4.5  Perveance 


Perveance  is  defined  as 


where  I  is  the  total  beam  current,  and  VCA  the  cathode  to  anode  voltage. 

The  total  beam  current  I  is  the  summation  of  current  emitted  from 
each  individual  ray  and  is  given  by 


I  -  l  J  CAREA 


where  JB  is  the  current  density  for  the  mth  diode,  and  CAREA^  Is  the 
ares  of  the  emitting  surface  of  the  mth  segment. 

Por  problems  represented  in  the  x  -  y  plane,  CAREAjj  is  just  the 
length  of  the  emitting  segment,  as  shown  in  Fig.  4.  For  problems  in  the 
z  -  r  plane,  CAREA^  is  a  figure  of  revolution  about  the  z-axis,  as  shown 
in  Fig.  4b,  and  is  computed  by  the  equation  for  a  frustum  of  a  right 
circular  cone. 


CAREA  -  w(Rl  +  R2) 

m  '  ' 


(a  +  b) 
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where  R1  le  the  radius  of  the  lower  base,  and  R2  is  the  radius  of  the 
upper  base*  This  foraula  was  chosen  so  that  when  deteraining  the  area 
of  a  spherical  or  linear  cathode,  the  same  foraula  could  be  used,  thus 
keeping  the  prograaalng  logic  sinple. 
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V.  RAY  TRACING  AND  SPACE-CHARGE  ASSIGNMENT 


Once  the  progrea  has  solved  Poisson's  equation  and  has  coaputed 
current,  its  final  task  is  to: 

1.  Coapute  the  particle  positions  of  the  electron  froa  the 


Lorentz  force  equation. 

2.  Assign  charge  to  the  aesh. 


5.1  Relativistic  Equation  of  Motion 


The  equations  of  aotlon  are  derived  froa  the  Lorentz  force  equa¬ 


tion. 


-e[E+M] 


where  e  is  the  charge  of  the  electron,  E  is  the  electric  field  obtained 

♦ 

froa  the  solution  of  Poisson's  equstlon,  v  is  the  velocity  of  the  elec- 
♦ 

tron,  end  B  is  the  external  aagnetlc  field  given  by  the  user.  The 
particle's  relativistic  aechanlcal  aoaentua  a ay  be  expressed 


♦  ♦  ♦ 

P  •  au  ■  ayv 


where  u  is  a  pseudo  particle  velocity  and  y  is  the  relativistic  factor 


-B)' 
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Separating  Eq.  54  Into  its  coaponents  in  cylindrical  coordinates,  ve 

have 


.  2 
•  u  Tv 

It1  “  ‘  n[£r  +  veB*]+  r  (57a) 


3u 

5T  *  "  "tE,  -  Vrl  (57b> 


where  n  ■  e/a  (the  derivation  of  Eq.  57  is  given  in  Appendix  A).  The 
last  tera  of  Eq.  57a  represents  the  centrifugal  acceleration,  which 
results  froa  the  notion  of  the  particle  about  the  axis.  Note  we  have 
aasuaed  that  the  problea  has  cylindrical  symmetry  and  the  magnetic  field 
has  only  radial  and  axial  coaponents. 

The  9  component  of  the  velocity  is  computed  by 


v 


9 


-  r9 


B 


zc 


(58) 


where  Bzc  and  r£  are  the  axial  magnetic  field  and  radial  position, 
respectively,  at  the  starting  position  of  the  trajectory.  Equation  58 
is  known  as  Busch's  theorem,23  which  makes  it  possible  to  determine  the 
angular  velocity  of  an  electron  about  the  axis  of  a  symmetric  system 
without  requiring  a  detailed  knowledge  of  the  form  of  the  electric 
fields. 


5.2  Electric  Fields 

The  electric  fields,  used  in  the  Lorentz  force  equation,  are 
computed  froa 
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The  coaponents  of  the  electric  field  in  cylindrical  coordinates  assua- 
lng  9  syne  try  are 


Er  3r 


_  3V 

E  ■  -  r— 
z  3z 


To  coapute  the  electric  field  at  point  P  of  Fig.  5,  the  first  step  is  to 


calculate  the  first  and  second  derivatives  at  the  nearest  node.  In  Fig. 


5,  the  nearest  node  is  labeled  0,  and  the  six  adjacent  nodes  are  labeled 


1  through  6.  Expanding  in  a  Taylor  series  about  node  0  to  each  of  the 


six  neighbors,  the  following  set  of  equations  is  generated: 


Fig.  5. 


Point  P  in  relation  to  its  nearest 
node  0  and  its  six  neighbors. 
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This  system  of  equations  is  overdetermined,  since  there  are  five  unknown 


derivatives  and  six  equations.  A  method  for  the  treatment  of  overdeter¬ 


mined  linear  equation  systems  is  the  method  of  least  squares. 


A  residual  vector  r  may  be  defined  by  rearranging  Eq.  61, 
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which  we  would  wish  to  be  zero.  The  result,  in  general,  will  not  be 


zero,  and  the  least  square  solution  of  Eq.  61  minimizes  the  Euclidean 


norm  of  the  residual  vector. 


After  the  derivatives  at  node  zero  have  been  computed,  the  components  of 


the  electric  field  at  point  P  are  calculated  by 
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5.3  Magnetic  Fields 


The  program  accepts  values  of  B  measured  experimentally  along  the 


axis  from  which  off  axis  values  of  Bz  and  Br  are  calculated.  We  can 


express  the  magnetic  flux  density  in  terms  of  a  magnetic  vector  poten¬ 


tial  , 


B  ■  V  *  A 


and  expanding  in  axisymmetric  cylindrical  coordinates,  we  have 
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and 


1  a(rA9) 

z  “  r  3r 


(67b) 


The  magnetic  vector  potential  A0(r,  z)  can  be  computed  from 


AQ(r,  z)  -  -  /  rBz(r,  z)  dr 


(68) 


Expanding  the  magnetic  flux  density  in  a  Taylor  series  near  the  axis  and 
assuming  the  second  term  is  negligible  with  respect  to  the  first  term, 
integrating  Eq.  68  results  in 


AQ(r,  *) 


rB  (0,  z) 
z 

2 


(69) 


Substituting  Eq.  69  into  Eq.  67,  the  components  of  the  magnetic  flux 
density  are 


B 

r 


3B  (0,  z) 
r  z 

2  3z 


(70) 


and 


Bz  -  B(0,  z) 


(71) 


where 


B(0,  z) 


+  (z 


'i+1 


'i+1 


(72) 
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The  time  integrator  used  in  this  program  is  a  finite-difference 
approximation  leapfrog  scheme.24  This  method  leapfrogs  positions  and 
velocities  in  time,  as  shown  in  Fig.  6.  Positions  and  fields  are 
defined  at  integral  time  levels  (t  *  0,  DT,  2DT,  3DT,  ...)  and  veloci¬ 
ties  are  defined  at  half-integral  time  levels  (t  *  1/2DT,  3/2DT,  ...). 

Positions  and  velocities  are  obtained  by  integrating 


dr 

dt 


dz 

dt 


(74a) 


and 


dv 


dv 


dt 


m 


dt 


(74b) 


where  Ff  and  Fz  are  the  forces  for  an  electron  with  charge  -e  moving  in 

♦  *► 
an  electric  field  E  and  a  magnetic  field  B. 
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Fig.  6.  The  equations  of  motion  are  integrated  forward  in  time 
using  the  second  order  leapfrog  scheme .  Positions  at 
time  level  n  -  1  are  updated  using  velocities  at  time 
level  n  -  1/2,  velocities  at  time  level  n  -  1/2  are 
updated  using  forces  at  time  level  n,  and  so  forth. 


yvem 


7  -  -  e(E  +  vaB  )  +  — 

r  v  r  6  zJ  r 


F  -  -  e(E  -  v  B  ) 
z  '  z  0  r' 


The  pseudo  velocities  are  predicted  at  the  (n  +  l)th  time  step  by 


(75a) 


(75b) 


n+1/2 


n+1/2 


fn  n+1^ 
n-1/2  lFr  +  Fr  J 
Ur  +  2^ -  At 

(Fn+Fn+1) 

u"-l/2  +  L*_.  «  . 


(76a) 


(76b) 


The  prograa  Chen  averages  Che  pseudo  velocities  and  computes  the  true 
velocity  of  Che  particle. 


vn 

r 


1 

7 


-n+l/2  t  un-l/2-| 


(77a) 


where 


-n+l/2  .  n-1/2 

11  +  It 


(77b) 


Y n  ■  1  +  !L_  v(rn,  z°) 

c 


and 


-n+l/2  n-1/2  FnAt 

u  *  u  +  — 

m 


(78) 


(79) 


v(rn,  zn)  in  Eq.  78  is  the  potential  at  the  electron  position  (rn,  zn) 
and  u  of  Eq.  79  is  the  first  order  approximation  of  the  pseudo  velocity 
which  is  used  in  Eq.  77. 


The  electron  positions  are  now  updated  by 


This  simple  second  order  scheme  was  chosen  because  it  is  a  good 
compromise  between  accuracy,  stability,  and  efficiency.  The  compromise 
between  accuracy  and  efficiency  can  be  altered  in  two  ways:  either  by 
using  a  high  order  scheme  and  larger  time  step,  or  by  using  a  low  order 
scheme  and  smaller  time  step.  The  higher  order  schemes  need  force 
values  at  several  time  steps,  thus  having  more  operations  per  time 
step.  The  leapfrog  scheme  is  simple,  and  is  more  efficient  and  has  less 
storage  limitations  than  a  higher  order  method.  The  accuracy  of  this 
scheme  is  more  than  sufficient  to  integrate  the  Lorentz  force  equation. 


5.5  Initial  Conditions  for  the  Equations  of  Motion 

Initial  conditions  for  the  electron  trajectories  are  determined 
from  space-charge  limited  parallel  plane  diode  theory. 

The  starting  position  of  the  trajectory  is  the  point  halfway 
between  nodes  NAC  and  NAC-1  and  is  denoted  by  (rQ,  zq)  in  Fig.  7.  It  is 
assumed  that  the  electrons  initially  proceed  normal  to  the  cathode  until 
they  reach  (rQ,  zc)  •  From  Child's  law  and  the  kinetic  energy  equation, 
the  time  needed  to  reach  this  point  is 


t 

o 


3d 


1/3 


(2nC) 


^1/2 


(81) 


where 


d 


z  NAC 


(82) 
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The  tiae  step,  At,  used  In  integrating  the  equation  of  motion  ia  calcu¬ 
lated  by 


( t  frdt) 

o _ _ 

(NAC  -  1) 


(83) 


where  FRDT  is  a  multiplicative  constant  by  which  At  may  be  fractionally 
adjusted. 


The  initial  position  (r^,  z^)  and  initial  velocity  (uro>  UZQ) 


are,  respectively,  calculated  by 


z  »  z 


r  +  t"  |  ■=■  nC|  sin  0 
c  o 


-  3t 


zo  o 


ro 


o  +  ‘o  [f  "6] 

[H 
M 

3to  [I 


3/2 


cos  0 


3/2 


3/2 


cos  6 


3/2 


sin  0 


(84a) 


(84b) 


(84c) 


(84d) 


where  zc  and  r£  are  the  coordinates  of  the  emitting  points,  and  0  is  the 
angle  of  the  line  connecting  the  nodes  1  through  NAC. 


5.5.1  Choice  of  Time  Steps 

The  program  has  the  capability  to  change  the  time  step,  At,  when 
the  magnetic  field  becomes  strong  enough.  The  choice  of  the  time  step 
must  be  related  to  the  characteristic  physical  frequency  of  the  problem, 
which  in  this  case  is  the  Larmor  period, 
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(85) 


c  m 


To  accurately  reaolve  a  La r»o r  period,  there  should  be  at  least  30  time 
steps  per  period.  When  a  magnetic  field  Is  present,  the  code  will  check 
to  see  If  there  are  at  least  40  time  steps  per  period.  If  there  are  not 
enough  time  steps,  At  will  be  recomputed  as 


new  2* 

1  40  nB 


and  the  unit  charge,  DQ,  of  Eq.  88  Is  recomputed  as 


_^new  __old 
DQ  "  DQ 


“ .  new “ 
At _ 

.  old 

J 


5.6  Charge  Assignment 


The  charge  distribution  of  particles  whose  positions  vary  continu¬ 
ously  with  time  is  replaced  by  a  finite  set  of  charge  density  values. 
After  each  time  step  in  the  integration  of  the  equations  of  motion,  a 
unit  of  charge,  DQ,  is  assigned  to  the  3  nodes  of  the  triangle  which 
contains  the  electron  at  that  time.  This  sharing  of  charge  with  the 
three  nodes  "smooths"  the  finite  charge  density  distribution  and  weakens 
the  mesh  dependence  of  the  problem.25 

The  unit  of  charge,  DQ,  which  is  assigned  to  the  mesh  in  the  x  -  y 
plane  is  computed  from 
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(88a) 


-I  At 

(qstep  co) 


and  in  the  z  -  r  plane, 


00  "  [QSTEP  2x  eQj  (88b) 

where  1  is  the  individual  current  of  the  ray,  At  is  the  time  step,  e  is 

o 

the  permittivity  of  free  space,  and  QSTEP  is  a  program  parameter  which 
will  be  discussed  later.  For  trajectories  which  proceed  along  the  z- 
axis,  radial  weighting  is  removed  and  DQ  is  calculated  as  if  it  were  an 
x  -  y  plane  problem. 


“>  -  iQstih;;  (88c> 

The  charge  DQ  is  assigned  to  the  nodes  of  the  mesh  by  the  following 
method.  Each  trajectory  segment  of  length  A  is  divided  into  QSTEP 
segments  in  which  DQ  is  positioned  at  the  center  of  each  minor  segment, 
as  shown  in  Fig.  8.  After  each  minor  time  step  in  the  trajectory  calcu¬ 
lation,  DQ  will  be  in  a  triangle  of  area  Atotaj.  The  nearest  node 
routine  determines  which  triangle  contains  the  electron  (this  routine  is 
described  in  detail  in  True's  thesis26).  Using  the  electron's  position 
and  the  3  corners  of  the  triangle's  vertices,  the  triangle  is  subdivided 
into  3  triangles  of  area  Aj  ,  A2 ,  A^ ,  as  shown  in  Fig.  9.  Introducing 
the  charge  assignment  function,  V^,  the  source  term,  Sit  associated  with 
these  3  nodes  will  have  DQ  added  to  them  by  the  amount 
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where 


W.  -  -r — - —  (90) 
1  total 

and  i  ■  1,  2,  3.  Note  that  as  the  electron  gets  closer  to  a  particular 
node,  more  charge  will  be  assigned  to  that  mesh  point,  since  the  charge 
assignment  function  is  larger  for  that  node. 

Once  all  the  charge  has  been  assigned  for  each  electron  trajec¬ 
tory,  nodes  along  the  z-axis  will  have  their  radially  weighted  charge 
removed  by 


C  a 

B  C 

r  6 

I 

i-1 

z-r 

(z-axis) 

(z-axis) 

*  6  * 
l 

i-1 

x-y 

(91) 


where  w^  are  the  coupling  coefficients  determined  by  Eq.  15.  Thus,  the 
charge  laid  down  on  the  axis  is  treated  as  if  it  were  an  x  -  y  plane 
problem. 
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VI.  EVALUATION  OF  THE  ELECTRON  TRAJECTORY  PROGRAM 


In  this  section,  the  validity  of  the  deformable  mesh  analysis  to 
space-charge  flow  problems  is  investigated.  Three  space-charge-limited 
electron  flow  problems  are  considered  and  the  results  of  the  program  are 
compared  with  theoretical  and  experimental  results. 


6.1  Pierce  Gun 


The  first  space-charge  flow  problem  considered  is  the  rectilinear- 
flow  Pierce  gun,  shown  in  Fig.  10.  This  case  is  ideally  suited  in  test¬ 
ing  the  electron  gun  program,  since  the  computed  values  of  voltage, 
space-charge,  electron  trajectories,  current  density,  and  perveance  can 
be  compared  to  theoretical  values. 

In  a  Pierce  gun,  the  planar  diode  is  truncated  and  two  focus 
electrodes  are  positioned  along  the  beam's  edge  in  order  to  establish 
the  four-thirds  power  variation  in  voltage  with  distance.  The  voltage 
distribution  in  the  charge-free  region  is  given  by27 


■  v<{!)  cos  (5 9) 


(92) 


where  V  ,  R,  d,  and  6  are  defined  in  Fig.  10.  The  voltage  inside  the 
beam  is  determined  from  space-charge-limited  parallel-plane  diode  theory 
as  is  given  by 


JLl 


mi 


>V*», 


...  .y.-* 

»  *J«  *  •» *  »  • 


•  *  *.•  -  1 

'  \  ■  w*  «*  ’ 

v  V  *  ** 

«  m  "  a 


■  V 

V 

• 


For  the  Pierce  gun  solved  by  computer  analysis  in  the  x  -  y  plane,  the 
cathode  and  focus  electrode  are  fixed  at  zero  volts,  the  anode  is  held 
at  1000  volts,  and  the  cathode  to  anode  spacing  is  1  centimeter.  The 
cathode  is  modeled  by  12  elementary  parallel-plane  diode  segments  and 
the  near-cathode  region  extends  three  nodes  in  front  of  the  cathode 
surface. 

The  logical  diagram  of  the  Pierce  gun  generated  by  AUTOMESH  is 
shown  in  Fig.  11  and  the  resulting  relaxed  mesh  is  shown  in  Fig.  12. 
The  results  of  the  deformable  mesh  solution  are  given  in  Table  1,  and  a 
plot  of  electron  trajectories  and  equipotentials  are  shown  in  Fig.  13. 
The  values  of  perveance,  average  current  density,  and  radial  deviation 
of  electron  trajectories  from  the  computer  analysis  agree  well  with 
theoretical  values.  The  values  of  current  density  calculated  for  each 
emitting  node  are  in  close  agreement  with  those  predicted  from  theory, 
except  for  the  upper  edge  node.  The  cause  of  the  errors  associated  with 
the  upper  emitting  node  results  from  the  discrete  computer  modeling 
along  the  beam's  edge  boundary.28 

For  the  edge  trajectory  of  the  Pierce  gun,  the  space-charge  is 
assigned  to  nodes  along  the  beam's  edge.  Correct  space-charge  assign¬ 
ment  for  the  beam  edge  would  assign  charge  to  subcells  4,  5,  and  6  and 
no  charge  to  subcells  1,  2,  and  3,  as  shown  in  Fig.  14a.  However,  the 
program  assigns  a  single  value  of  space-charge  to  the  node,  which 
results  in  the  space-charge  density  being  uniformly  distributed  through¬ 
out  the  dodecagon  area  surrounding  the  node,  as  shown  in  Fig.  14b.  As  a 
result,  the  space-charge  density  is  lower  on  the  beam  edge,  which  causes 
the  potentials  to  be  higher  than  those  within  the  beam.  Since  the 
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Table.  1.  Computer  results  for  the  planar  Pierce  gun 


Theo reties 1 

Calculated 

Percent 

Value 

Value 

Error 

Current  Density 

7. Ml  x  10"21 

7.435  x  10'2 

0.74X 

Perveance 

2.334  x  10'8 

2.351  x  10‘8 

0.74X 

Current  Density  At  The  Emitting  Nodes 


Emitting 

Node 

Theoretical  Value 
(anp/  sq.  m) 

Calculated  Value 
(amp /  sq.  m) 

Percent 

Error 

1 

7.381  x  10‘2 

7.52  x  10‘2 

1.889 

2 

7. Ml  x  10'2 

7.436  c  10'2 

0.753 

3 

7. Ml  x  10‘2 

7.378  x  10‘2 

-0.033 

4 

7. Ml  x  10'2 

7.402  x  10'2 

0.289 

5 

7. Ml  x  10‘2 

7.413  x  10‘2 

0.435 

6 

7.381  x  10‘2 

7.422  x  10‘2 

0.567 

7 

7. Ml  x  10'2 

7.436  x  10*2 

0.751 

8 

7.381  x  10*2 

7.446  x  10'2 

0.681 

9 

7.381  x  10'2 

7.476  x  10‘2 

1.294 

10 

7. Ml  x  10*2 

7.299  x  10‘2 

-1.109 

11 

7.381  x  10'2 

7.298  x  10'2 

-1.119 

12 

7. Ml  x  Kf2 

8.047  x  10'2 

9.02 

Maximum  Radial  Deviation  In  Electron  Trajectories 


Emitting 

Node 

Theoretical 
Deviation  (cm) 

Calculated 
Deviation  (cm) 

Percent 

Error 

1 

0.0 

0.0 

0.0 

2 

0.0 

-4.20  x  10'4 

-0.464 

3 

0.0 

1.2  x  10’3 

0.710 

4 

0.0 

3.727  x  10*4 

0.137 

5 

0.0 

8.636  x  10*4 

0.237 

6 

0.0 

5.545  x  10"4 

0.122 

7 

0.0 

1.045  x  10° 

0.191 

8 

0.0 

8.364  x  10*3 

0.131 

9 

0.0 

2.327  x  10*3 

0.319 

10 

0.0 

-2.182  x  10'3 

-0.267 

11 

0.0 

1.009  x  10'3 

0.122 

12 

0.0 

2.727  x  10'5 

0.003 
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Fig.  15.  Relaxed  mesh  with  a  beam-edge  form- 
line  for  the  planar  Pierce  gun. 
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Table  2.  Computer  results  for  the  planar  Pierce 
gun  using  the  beam  edge  guard  mesh. 


Perveance 


Theoretical 

Calculated 

Value 

Value 

7.381  x  10‘2 

7.409  x  10"2 

2.334  x  10-6 

2.343  x  10"6 

Percent 

Error 


Emi tti ng 
Node 


Theoretical  Value 
(amp  /  sq.  m) 

Calculated  Value 
(amp/  sq.  m) 

7.381  x  10”2 

7.520  x  10‘2 

7.381  x  10"2 

7.438  x  10'2 

7.381  x  10"2 

7.382  x  10"2 

7.381  x  10"2 

7.406  x  10"2 

7.381  x  10"2 

7.416  x  10'2 

7.381  x  10"2 

7.427  x  10~2 

7.381  x  10"2 

7.438  x  10"2 

7.381  x  10'2 

7.451  x  10'2 

7.381  x  10"2 

7.477  x  10'2 

7.381  x  10‘2 

7.348  x  10‘2 

7.381  x  10'2 

7.2  x  10"2 

7.381  x  10'2 

7.52  x  10"2 

1.86% 

0.78% 

0.01% 

0.34% 

0.48% 

0.63% 

0.78% 

0.95% 

1.3% 

0.44% 

2.45% 

1.88% 
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Another  case,  which  is  also  Ideally  suited  in  testing  the  electron 
gun  progress  is  the  spherical  diode.29  The  spherical  diode  consists  of 
two  concentric  spheres  in  which  the  inner  surface  of  the  outer  sphere  is 
coated  with  an  eaisslve  material  and  the  inner  sphere  is  made  positive. 

The  current,  which  by  symmetry  is  everywhere  radial,  flows  between 
the  spheres  and  is  given  by 


29.34  x  1Q~6  y3/2 


total 


where 


2  3  4 

a  -  y  -  0.3  y  +  0.075  y  -  0.00143  y  +  ... 


fe) 


y  -  In 


Figure  17  shows  the  geometry  of  the  diode,  which  consists  of  a  Neumann 
boundary  along  the  z-axis.  Since  the  problem  is  axisymmetric,  the  two 
half-circle  boundaries  are  actually  figures  of  revolution  about  the  z- 
axis.  The  logical  space  diagram  for  this  example  is  shown  in  Fig.  18, 
and  the  relaxed  mesh  is  shown  in  Fig.  19. 
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Table  3  compares  the  results  of  the  computer  analysis  with  that  of 
the  theoretical  analysis,  and  the  plot  of  the  electron  trajectories  and 
equipotentials  are  shown  in  Fig.  20.  Again,  the  values  of  emitted 
current  and  radial  deviation  of  electron  trajectories  produced  by  the 
analysis  of  the  computer  agree  very  well  with  that  of  theory. 


Table  3,  Computer  results  for  the  spherical  diode. 
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Theoretical  Calculated  Percent 

Value  Value  Error 
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Perveance  (microperus) 


Average  Error  in  Current 
Density  At  The  Cathode 
Surface  (amps/sq  m) 

Average  Error  In  Radial 
Deviation  of 
Trajectories 


11.675 


11.7366 


0.53% 


0.91% 


1.18% 
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Equlpotentlals  and  electron  trajectories  for  the  spherical  diode 


6.3  Mod-I  Gun 


In  each  of  the  cases  discussed  above,  analytical  solutions  were 
available  to  verify  the  results  of  the  computer  program.  However,  for 
guns  of  practical  interest,  analytic  solutions  are  not  available. 
Therefore,  a  comparison  of  computed  and  experimental  results  for  actual 
guns  are  of  interest.  In  this  section,  evaluation  of  the  simulation 
capabilities  of  the  computer  program  is  conducted  using  an  experimental 
gun  design  (MOD-I)  developed  at  Bell  Telephone  Laboratories  by  R.  D. 
Brooks  (November  1969). 

The  dimensions  of  the  MOD-I  gun  design  are  shown  in  Fig.  21.  The 
geometry  of  this  gun  is  shown  in  Fig.  22,  and  the  logical  diagram  is 
shown  in  Fig.  23.  The  resulting  relaxed  mesh  is  given  in  Fig.  24,  and 
Fig.  25  shows  the  plot  of  equipotentials  and  electron  trajectories. 
Note,  the  cathode  is  graded  such  that  the  diode  segment  lengths  are 
smaller  towards  the  end  of  the  cathode.  This  grading  will  better  model 
the  cathode,  since  a  large  percentage  of  the  total  beam  current  is 
emitted  from  the  cathode  edge. 
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Comparisons  between  calculated  and  experimental  beam  parameters 


are  tabulated  in  Table  4.  The  accuracy  of  the  location  of  the  beam 
minimum  and  the  value  of  the  perveance  computed  by  the  program  is  in 
good  agreement  with  the  experimental  results.  The  discrepancy  between 
the  experimental  and  computer  results  might  be  in  part  due  to  experimen¬ 
tal  measurement  error  and  small  dimensional  differences  between  the  gun 
tested  by  Brooks  and  the  gun  dimensions  used  in  the  computer  analyses. 
The  discrepancy  in  the  measured  and  calculated  beam  diameter  at  the  beam 
minimum  is  probably  due  to  thermal  effects.  The  problem  of  accurately 
computing  the  effects  of  thermal-velocity  spread  at  the  cathode  is  not 
modeled  in  the  computer  program. 


Table  4.  Summary  of  results  for  the  2D  axisymmetric  MOD-I  gun. 


Experimentally 

Measured 

Computer 

Results 

Percent 

Error 

Perveance 

(microperv) 

1.63 

1.604 

1.6% 

Beam 

Minimum 

Location 

(cm) 

6.655 

6.833 

2.5% 

Beam 

Minimum 

Diameter 

(cm) 

0.8585 

0.9305 

8.4% 

VII.  CONCLUSIONS 


As  was  shown  in  the  preceding  section,  agreement  with  theoretical 
and  experimental  results  with  those  of  the  computer  program  is  very 
good.  Perveance  error  is  in  the  one  percent  range,  and  error  in  current 
density  along  the  emitting  surface  is  typically  a  few  percent  or  less. 
The  results  of  Chapter  VI  indicate  that  the  overall  performance  of  the 
deformable  mesh  program  is  very  good. 

Some  areas  for  further  work  to  improve  the  computer  program  are: 

1.  Accurate  representation  of  magnetic  fields.  The  accurate 
representation  of  such  fields  might  be  improved  by  including 
higher  order  terms  in  the  Taylor  series  expansion,  or  using  an 
elliptic  integral  method  for  representing  axisymmetric  mag¬ 
netic  fields.30  Another  alternative  would  be  to  use  the 
POISSON  group  programs  to  solve  the  magnetic  fields.  These 
sets  of  programs  were  originally  written  to  solve  the  magneto¬ 
static  vector  potential  problem  using  nonlinear  iron. 

2.  As  mentioned  in  Section  6.3,  the  effect  of  thermal-velocity 
spread  at  the  cathode  is  not  included  in  the  program.  A 
method  of  modeling  thermals  could  be  incorporated  into  the 
code. 

3.  Since  determining  the  positions  of  grid  wires  is  tedious,  the 
code  could  be  modified  to  generate  spherical  grid  data,  addi¬ 
tional  guard  lines,  and  computation  of  grid  interception. 
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APPENDIX  A 


DERIVATION  OF  THE  EQUATIONS  OF  MOTION 

The  position  of  a  particle  can  be  described  in  cylindrical  coordi¬ 
nates  R,  0,  a  by  the  position  vector 

♦  *  * 

r  "  Rr  +  zz  (A.l) 

A  A 

where  r  is  a  unit  vector  in  the  x  -  y  plane  and  z  is  the  unit  vector  in 
the  z  direction. 

The  velocity  and  acceleration  vectors  are  found  by  differentiating 
Eq .  A.l, 

A  A 

V-dr»Rr+r^~+zz+z^  (A.  2) 

Note,  this  invovles  derivatives  of  unit  vectors  which  are 

A 

“  08  (A. 3a) 


and 


de 

dt 


(A. 3b) 


The  unit  vector  z  does  not  change  in  direction,  so  its  time  derivative 
ia  zero.  Substituting  Eq.  A. 3a  into  Eq.  A. 2,  the  velocity  vector  is 
given  by 


(A. 4) 


V  -  Rr  +  V.0  +  zz 

where 


V  - 

o 


R0 


(A.  5) 


Differentiating  Eq.  A. 4  results  in  the  acceleration  vector 


+ 

a 


[R  -  r02]r  +  [R0  +  V  ]  + 


•  •  A 
22 


(A. 6) 


The  equations  of  motion  are  derived  from  the  Lorent2  force  equation, 


dP 

dt 


V  X 


B] 


(A. 7) 


where  the  relativistic  mechanical  momentum  is  given  by 


>  ♦  ♦ 
P  «  uryv  =  mu 


(A. 8) 


The  velocity  vector  u  is  a  pseudo  particle  velocity  where 


dr 

dt 


y  v 


(A. 9) 


where  y  is  the  relativistic  factor  given  by 


y  m 


(A. 10) 


Substituting  Eq.  A. 8  into  Eq.  A. 7  gives 
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(A. 11) 


♦  r  +1 

-  -  —  E  +  JLJL1 
it  m  T  I 


Setting  Eq.  A. 6  equal  to  Eq.  A. 11  and  separating  into  its  cylindrical 
components  results  in  the  relativistic  equations  of  motion,  which  are 
used  in  the  computer  program  to  solve  for  the  particle's  position  and 


velocity. 


uB  "1 
0  z  •  . 

-  n  Er  +  — —  +  r0 

_ 


r  uaB 

.  „  .  0  r 

n  E  +  - 

z  Y 


du„  TuB  u  B  1  uu„ 

_ 9  _  z  r  _  r  z  _  r  6 

dt  "  n  y  Y  Yr 


(A. 12a) 


(A. 12b) 


(A. 12c) 


where  n  =  e/m. 
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MISSION 

of 

Rome  Air  Development  Center 

RAVC  plan s  and  execute 6  tie.6za.Ach,  development,  tz6t 
and  ielected  acquisition  program 6  in  Support  c f 
Command,  Control,  Communications  and  Intelligence 
(Co  I)  activities.  Technical  and  engineering 
support  within  areas  of  competence  is  provided  to 
ESV  Program  Quizes  (PCs)  and  other  ESV  elements 
to  perform  elective  acquisition  c  <  C3 1  systems. 

The  areas  o &  technical  competence  include' 
communications ,  command  and  control,  battle 
management,  information  processing,  surveillance 
sensors,  intelligence  data  collection  and  handling, 
solid  state  sciences,  electromagnet ic s ,  and 
propagation,  and  electronic,  maintainability, 
and  compatibility. 


